clc; clear;
cd 'data'

lambda_data = csvread('lambda_m.csv');
rho_data = csvread('rho.csv');
Y_data = csvread('Y_m.csv');
beta = csvread('beta.csv');
rhodot = csvread('rhodot.csv');
sm = 1-rhodot;
 
cd ..

I = 44;

alpha = 1; 
theta = 2.5; sigma = 3;
YEAR = 15;

theta_n = 2.7; sigma_n = 2.8;  

lambda = reshape(lambda_data,I,I);    
rho = 1-reshape(rho_data(:,YEAR),I,I);   
rho_autarky =  max(0,1-(1-diag(rho))./diag(lambda));
 
Y = Y_data;
 
Gini = alpha.*(sigma-1)./(sigma*theta) - 2.*alpha.*(theta-(sigma-1))./(sigma*theta).*(sigma-1)./(4*theta-2*(sigma-1))...
    .* sum(lambda.*(1-rho).*repmat(Y',I,1)./repmat(Y,1,I),2);
Gini_autarky = alpha.*(sigma-1)./(sigma*theta) -  2.*alpha.*(theta-(sigma-1))./(sigma*theta).*(sigma-1)./(4*theta-2*(sigma-1)).* (1- rho_autarky);
 
Gini_M = beta.*(sigma-1)./(sigma*theta) + (1-beta).*(sigma_n-1)./(sigma_n*theta_n) - ...
(theta-(sigma-1))./(sigma*theta).*(sigma-1)./(2*theta-(sigma-1)).*(beta)./(sm).* sum(lambda.*(1-rho).*repmat(Y',I,1)./repmat(Y,1,I),2);

Gini_M_autarky = beta.*(sigma-1)./(sigma*theta) + (1-beta).*(sigma_n-1)./(sigma_n*theta_n) - ...
(theta-(sigma-1))./(sigma*theta).*(sigma-1)./(2*theta-(sigma-1)).*(beta)./(sm).*(1- rho_autarky);

K = 500;
grid = linspace(10^(-3),1-10^(-3),K);

for p = 1:K
   s = grid(p);
   DOM = 1./(1-diag(rho)).*(((1-s)./(1-diag(rho))).^((1-sigma)/theta) - 1);
   DOM = (s>=diag(rho)).*DOM;
   r(:,p) = (1-beta.*(sigma-1)./(sigma.*theta)-(1-beta).*(sigma_n-1)./(sigma_n.*theta_n)) + beta./sm.*sum((theta-(sigma-1))./(sigma.*theta).*DOM,2);
   LAMBDA = lambda.*repmat(Y',I,1)./repmat(Y,1,I)./(1-rho).*(((1-s)./(1-rho)).^((1-sigma)/theta) - 1);
   LAMBDA = max(0,abs(1-eye(I,I)).*LAMBDA);
   INT(:,p) = beta./sm.*(theta-(sigma-1))./(sigma.*theta).*sum(LAMBDA,2);
   %----------------------------------------------------------------------%
   DOM = 1./(1-diag(rho)).*(((1-s)./(1-diag(rho))).^((1-sigma_n)/theta_n) - 1);
   DOM = (s>=diag(rho)).*DOM;
   yn(:,p) = (1-beta)./(1-sm).*sum((theta_n-(sigma_n-1))./(sigma_n.*theta_n).*DOM,2) + (1-beta.*(sigma-1)./(sigma.*theta)-(1-beta).*(sigma_n-1)./(sigma_n.*theta_n)) ;
   %----------------------------------------------------------------------%
   DOM = 1./(1-rho_autarky).*(((1-s)./(1-rho_autarky)).^((1-sigma)/theta) - 1);
   DOM = (s>=rho_autarky).*DOM;
   ra(:,p) = (1-beta.*(sigma-1)./(sigma.*theta)-(1-beta).*(sigma_n-1)./(sigma_n.*theta_n)) + beta./sm.*sum((theta-(sigma-1))./(sigma.*theta).*DOM,2);
end

ym = (r+INT);
ya = ra;
Nm = floor(100*sm); Nn = 100-Nm;

for i = 1:I
    yy = [repmat(ym(i,:),Nm(i),1); repmat(yn(i,:),Nn(i),1)];
    yy = reshape(yy,100*K,1);
    Gini_Total(i,:) =  sum(sum(abs(kron(yy',ones(100*K,1)) - kron(yy,ones(1,100*K)))))./(2*(100*K)^2*mean(yy));
    yy = [repmat(ya(i,:),Nm(i),1); repmat(yn(i,:),Nn(i),1)];
    yy = reshape(yy,100*K,1);
    Gini_Total_Autarky(i,:) =  sum(sum(abs(kron(yy',ones(100*K,1)) - kron(yy,ones(1,100*K)))))./(2*(100*K)^2*mean(yy));
end

A = sortrows([Gini_M_autarky./Gini_M Gini_Total_Autarky./Gini_Total],2);
bar(A,'DisplayName','A')

ylim([0.2 1]);

xticks([1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44]);
xticklabels({"SVN","SVK","LTU","EST","HUN","IRL","NOR","BEL","CZE","DEU","BGR","NLD","HRV","LVA","LUX","AUT","CAN","POL","IDN","SWE","ROU","DNK","PRT","FIN","TUR","MEX","TWN","RUS","MLT","CHE","ITA","KOR","ESP","FRA","GRC","BRA","GBR","IND","ROW","AUS","CHN","JPN","CYP","USA"}); 






